Tight-binding study of structure and vibrations of amorphous silicon 
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We present a tight-binding calculation that, for the first time, accurately describes the structural, 
vibrational and elastic properties of amorphous silicon. We compute the interatomic force constants 
and find an unphysical feature of the Stillinger- Weber empirical potential that correlates with a much 
noted error in the radial distribution function associated with that potential. We also find that the 
intrinsic first peak of the radial distribution function is asymmetric, contrary to usual assumptions 
made in the analysis of diffraction data. We use our results for the normal mode frequencies and 
polarization vectors to obtain the zero-point broadening effect on the radial distribution function, 
enabling us to directly compare theory and a high resolution x-ray diffraction experiment. 

PACS numbers: 61.43.Dq, 62.20.Dc, 63.50.+X, 78.55.Qr 
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Amorphous silicon (a-Si) is a prototype for continuous- 
random-network covalent glasses that, with some hydro- 
gen content, has technological applications as a relatively 
inexpensive electronic material. While the basic struc- 
ture of a-Si is believed to be a four-fold-coordinated con- 
tinuous random network, detailed information about net- 
work connectivity and defects is lacking. Atomic resolu- 
tion structure is very difficult to determine directly, and 
experiments have relied on unusual or indirect probes 
such as variance coherence microscopy Q and Raman 
spectroscopy 0, as well as on more standard tech- 
niques such as diffraction 

HE1 and EXAFS 0,0. The 
experimental measurements suggest significant deviation 
from a continuous random network, including average co 
ordination that is significantly less than 4 (e.g. Rcf 
and that unanncaled samples may be paracrystallinc 
Many empirical-potential simulations have been done, 
but it is not clear if empirical potentials are accurate 
enough to give reliable results for properties, such as co- 
ordination defects, that depend on bond breaking and 
bond formation. A number of simulations of a-Si struc- 
ture have used electronic-structure based methods, which 
are generally among the most reliable for solid state sys- 
tems (e.g. Refs. 0, la E3, E3) 

However, none have care- 
fully compared the radial distribution function (RDF) 
to high resolution experiments and none included 
quantum- mechanical vibrational effects. Another impor- 
tant question concerns the vibrational properties of a-Si, 
which give us information about the structure and the 
interactions of atoms in the material. The vibrational 
density of states (VDOS) was measured experimentally 
using inelastic neutron scattering (INS) Empirical- 
potential simulations have been used to analyze vibra- 
tional properties in detail but all show significant 
errors in the shape of the VDOS or in other properties. 
While the VDOS of a-Si has been simulated with elec- 
tronic structure methods H, 0, Q] , the underlying force 
constants themselves have not been analyzed. There have 
been many studies of force constants in crystalline Si, 
which shows unusual phonon dispersion and force con- 



stants that oscillate in magnitude as a function of dis- 
tance [HE3- 

We study the elastic constants, vibrational properties, 
and structure of a-Si using a tight-binding (TB) total- 
energy method. We find elastic constants and VDOS that 
are in good agreement with experiment, and qualitatively 
better than empirical-potential simulations. The struc- 
ture has a sharp first-neighbor RDF peak that agrees 
very well with experiment when zero-point and thermal 
broadening is included. This peak is significantly non- 
Gaussian, calling into question the coordination-statistics 
analysis of previous diffraction experiments. 

We use the Naval Research Laboratory (NRL) TB 
method 0,^1- The non-orthogonal sp 3 -basis TB model 
has been shown to accurately describe the elastic con- 
stants and phonon dispersion in crystalline Si and the 
electronic density of states for a highly defected amor- 
phous model . To generate the a-Si models we re- 
lax using TB-calculated forces a-Si models generated 
by other techniques. The NRL-TB model is used to 
calculate the energy of the structure and the atomic 
forces jilj. The conjugate-gradient method is applied 
to find mechanical-equilibrium positions at a fixed vol- 
ume, employing the criterion that components of atomic 
forces be less than 10 -3 cV/A. The relaxation procedure 
is carried out at several volumes to obtain results at zero 
pressure, but components of the stress tensor, generally 
of magnitude less than 0.8 GPa, remain. 

One model, which we denote TBI, is generated by re- 
laxing (using TB)a 216 atom perfect continuous-random- 
network model |l9j with periodic boundary conditions 
relaxed with a Keating interatomic potential [2(|. The 
TB-relaxed model is perfectly four-fold coordinated, with 
1.3% lower density than the crystal, compared to 1.7% 
lower density measured experimentally p|. The bond- 
angle distribution has a RMS deviation of 11° from the 
average value of 109.2°, in close agreement with relaxed 
ab initio calculation |l0l | and analysis of experiment Q . A 
second model, which we denote TB2, is generated by re- 
laxing a structure from a molecular-dynamics simulation 
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TABLE I: Selected elastic constants c, bulk modulus B and 
Young's modulus E (10 11 dyn/cm 2 ). The index i varies from 
1 to 3, and j from 4 to 6. 





TBI 


TB2 


Exp./FP 


SW (a) 


Cji 


16.31-16.45 


15.06-16.00 


13.8 l,, ^17(2) (c, 


11.94-13.11 




5.68-5.84 


5.26-5.56 


4.8 (6) , 4.5 (a) 


2.54-3.21 


P 


5.77 


5.06 




2.62 


C12 


4.77 


5.32 




6.69 


B 


8.73 


8.99 


5.9 (e) ,8.25 (/) 


8.52 


E 


14 (s) 


13 (s) 


12.4(3) (a) 


7(9) 



11.7(5)-13.4(5) (fe) 



(a) Ref. 24; (b) Ref. ,25.; (c) Ref. 2£; 
(d) Defined here as (cn-ci2)/2; (e) Ref. ^7 
(f) Ref. L£; (g) based on values of C12 and c p ; 
(h) Ref. 28. 



of the rapid quenching of liquid Si using the environment 
dependent interatomic potential ■ The TB2 structure 
is slightly more dense than TBI, but still about 0.5% less 
dense than the crystal. The energy is 28 meV/atom lower 
than the TBI energy, despite the presence of 6% 5-fold 
and 0.46% 3-fold coordinated atoms (corresponding to an 
average coordination of 4.05) [2^. The RMS bond-angle 
deviation is 12.5°, although the distribution has wide, 
non-Gaussian wings; excluding 2% of the bond-angles re- 
duces the RMS deviation to 10.4°. We also show some 
comparisons with results using the Stillinger- Weber (SW) 
interatomic potential j2^. The SW potential, which in- 
cludes radial and bond-angle terms, is one of the most of- 
ten used potentials for simulations of Si. We use a struc- 
ture (Ref. Eil Table II, model IV) generated by relaxing 
with SW the same starting structure as TBI. Finally, we 
note that while it is possible to use electronic structure 
methods to generate amorphous structures from proce- 
dures that are less dependent on the initial structure, 
it is very expensive computationally. The difficulty in 
fully annealing the structure seems to lead to a consis- 
tent overestimate of the width of the first-neighbor peak 
in the RDF HQ. 

The relaxed static lattice TB elastic constants Cjj were 
obtained by the method of homogeneous deformation. 
The TB results j3]| are compared in Table [I] with results 
of first-principles (FP) calculations, SW calculations, 
and several experiments on dense samples (a wider range 
of shear values are quoted in Table V of Ref. l30|) . Al- 
though there is some deviation between the two TB struc- 
tures it is small. While ultrasonic measurements of elastic 
properties are not available for a-Si, the Young's modulus 
E can be measured with a vibrating reed apparatus, and 
other elastic constants can be inferred from spectroscopic 
studies. Our TB results for both models are close to the 
experimental values, although our value of C44 is likely 
10-20% too large. The SW empirical potential results 
are significantly worse in comparison with experiment. 

The VDOS is calculated from a dynamical ma- 
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FIG. 1: Vibrational density of states of a-Si. A Gaussian 
broadening of FWHM — 20 cm -1 is employed. The experi- 
mental data are from Ref. |32T [. 

trix approach. The matrix elements ^ a p(i,j) = 
AF a (i) / Aup(j) are calculated using the TB forces with a 
central-finite difference approach that eliminates all odd- 
order anharmonic terms in the potential energy [3j| . The 
TB VDOS for structures TBI and TB2 are compared 
with SW results and INS [33] measurements in Fig. ^ 
For both structures the TB calculation yields the overall 
shape very well; it exactly describes the low frequency TA 
peak, gives a slightly too small frequency of the LA peak 
(300 cm -1 ) and about a 10% percent too high frequency 
of the high frequency TO peak. The TB results are a 
qualitative improvement over results based on the SW 
potential, as shown in the figure, and they are in good 
agreement with ab initio results for a 216 atom structural 
model Q. 

The range of the effective interactions in the solid can 
give us information about the physics of the interactions, 
and can guide the development of approximations such as 
empirical potentials. In Fig. [21 we plot all of the cartesian 
force constants between pairs of atoms with interatomic 
distances less than 10 A. The difference in range between 
the SW results and the TB results is easy to see: The SW 
interactions are large up to about 3.5 A, and go to exactly 
zero at twice the SW cutoff of 3.75 A. The TB interac- 
tions are already quite small at 2.8 A, but do not go to 
zero even at 10 A. This comparison of TB and SW leads 
to a view of interactions in the solid that is more sub- 
tle than the usual assumption that empirical potentials 
are short ranged and that the real interactions are long 
ranged: The SW potential interactions go to zero at a 
range that is too short, but at intermediate distances the 
interactions are too strong. We also note that the pre- 
ponderance of force constants as a function of interatomic 
distance give a clear envelope function that has an oscil- 
latory behavior which matches the RDF peak positions. 
This is qualitatively similar to the case of the crystal, 
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FIG. 2: Force constants between pairs of atoms in SW (top) 
and TBI (bottom) relaxed structures (dots). Superimposed 
are the corresponding J(r) functions (jagged lines) in arbi- 
trary units and the experimental, annealed-sample, results of 
Ref. for J(r) (smooth lines). The upper scatterplot in the 
TB panel is a magnification of the smaller magnitude force 
constants. 



even though the explanations for the oscillation in the 
crystal do not apply to the amorphous structure 0, 0] . 

The problem with the SW potential is a direct conse- 
quence of the form of the potential. In the amorphous 
there are pairs of atoms in the second-neighbor peak with 
distances smaller than the SW cutoff. It is clear from 
the TB force constants that the effective interactions 
for these pairs is qualitatively different from the first- 
neighbor interactions. However, in the SW simulation 
these second-neighbor pairs interact through terms that 
are meant to describe the interactions of first-neighbor 
atoms. In particular, the two-body contribution has 
strong negative curvature at these distances, and the 
three-body terms include contributions from triplets with 
a vertex angle that does not correspond to an atom with 
two sp 3 orbitals in bonding configurations. These two 
types of contributions lead to the unphysically large force 
constants in the SW results at this range of distances. 
The range of incorrect force constants also coincides with 
the shoulder in the SW RDF that is not observed in our 
TB results or in the experimental measurements pjij . 

The distribution of force constants gives us information 
about the types of effective interactions between bonded 
atoms. Under the first peak of the RDF the largest pos- 
itive cartesian force constants are twice the magnitude 
of the largest negative force constants for both SW and 
TB. This relation is consistent with an effective bond- 
stretching interaction for first-neighbors. We plot the 
results for the bond-strctching components in a plot as 
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FIG. 3: (a) First peak of the static RDF and TBI bond 
stretching force constants, (b) Broadened results correspond- 
ing to T=10K in comparison with experiment (annealed 
sample) . 



a function of r (Fig. Ot). The radial force constants de- 
crease with increasing r as one expects from a physi- 
cally reasonable first-neighbor bonding potential. Pairs 
with large (small) interatomic bond stretching force con- 
stants will have small (large) relative mean square dis- 
placements, so these results clearly have an impact on 
the nature of the broadening of the RDF. 

Very little attention has been given in the literature to 
the shape of the first peak in the RDF J(r) [35j . This 
peak has been measured very carefully at T = 10 K with 
x-ray diffraction, using high energy photons and high res- 
olution, i.e., large Q max , by Laaziri et al. [jj. They obtain 
a fit of their data to a Gaussian, with average coordi- 
nation of 3.88±.0I (3.79±.01) for the annealed (unan- 
nealed) sample. In Fig. |3K we plot the first peak of the 
static J(r) for models TBI and TB2, and the SW re- 
sults. The TB static peak is asymmetric, and its width 
is significantly larger than the static-disorder estimate by 
Laaziri et al. In order to compare directly with the ex- 
perimental J(r) it is necessary to properly take account 
of the zero-point and thermal broadening. The quan- 
tity measured by the x-ray experiment is, in the small- 
displacement limit, 



J(r) 



1 N 

- y 



¥ exp(-(r ij -ry/(2Ul j )), 



where IK = 



{(vij ■ Uij) 2 }. Thus we need the mean- 
squared relative displacements, for pairs of atoms, along 
the interatomic vector direction. We calculate them 
within the harmonic approximation at T — 10 K using 
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our computed vibrational modes. Since T = 10 K essen- 
tially corresponds to T — K for these considerations, 
what we obtain is the minimum measurable width for 
the first peak in the RDF of amorphous silicon. As seen 
in Fig. 12)3 the results are in agreement with experiment, 
aside from a small skewing of the theoretical function to 
large r. Although it has not been observed in a-Si, this 
type of asymmetry has been observed in EXAFS of amor- 
phous germanium [2^. Both the TBI and TB2 models, 
despite the very different originating structure and differ- 
ences in coordination defects, show nearly identical RDF 
first peaks. The good agreement with experiment of the 
broadened RDF suggests that our static peak width is 
correct, and that Laaziri et al. underestimate the static 
disorder contribution to the broadening. This may be 
caused by inaccuracy in the polycrystalline J(r) that is 
used to estimate the dynamic broadening. In the experi- 
ments a lower Q max (35 A -1 ) was used for the polycrystal 
than for the amorphous structure (55 A^ 1 ), although the 
former is expected to have a narrower first peak. Numer- 
ous other treatments using EXAFS or diffraction have 
not been considered here because they all use too low 
values of Qmax for obtaining reliable information on the 
first peak. The only other theoretical study of quantum 
effects in J(r) is by Herrero who used the SW po- 
tential but treated the quantum-effects on the nuclear 
vibrations exactly. Our results using the SW potential 
are presented in Fig. [3] The result for the amount of 
zero-point broadening is consistent with Herrero's work, 
although due to differing approximations a direct com- 
parison is not possible. We note that the Wooten model 
on which both the SW and TBI models are based yields 
a static J(r) (not shown) that is quite symmetric, and as 
broad as the experimental breadth. 

To conclude, we have shown that the NRL-TB method 
can reliably compute structural, vibrational, and elastic 
properties of a-Si. The results are nearly identical for 
two structural models, one with perfect four- fold coordi- 
nation and one with several atomic percent coordination 
defects. We have presented the first discussion of force 
constants in a-Si, which has revealed limitations of the 
most frequently used empirical potential for silicon. Our 
calculated elastic constants fall within the range of ex- 
perimental values for imperfect samples prepared under 
various conditions. We have also carefully studied the 
first peak in the radial distribution function. Wc observe 
a clear asymmetric peak in the case of the static quan- 
tity which is not observable experimentally We have 
included the (essentially) zero-point broadening effects 
in J(r) to obtain the experimentally measured quantity. 
Our two structural models, which have average coordi- 
nations of 4.00 and 4.05, respectively, reproduce the first 
peak in the experimental J(r) (for the annealed sam- 
ple) except for a slight asymmetry still present in the 
broadened result. We believe that such an asymmetry 
is expected on physical grounds and that perhaps it has 



been "missed" experimentally because of the challeng- 
ing analysis required to obtain J(r) from the diffraction 
results. 
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